% Copyright (C) 2009,2010,2011,2012  Marco Restelli
%
% This file is part of:
%   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
%
% LDGH is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% LDGH is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
% License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LDGH. If not, see <http://www.gnu.org/licenses/>.
%
% author: Marco Restelli                   <marco.restelli@gmail.com>


function ui = dgcomp_interpolation( grid1,ddc_grid1,base1,uuu1, ...
                                    xxx )
% ui = dgcomp_interpolation( grid1,ddc_grid1,base1,uuu1, ...
%                            xxx )
%
% Interpolate a DG-NS solution from grid1 to the unstructured points
% xxx (as column vectors).
%
% The input arguments grid1, ddc_grid1 and uuu1 can be arrays in order
% to deal with domain decomposition (base1 is a scalar since we assume
% that all the partitions in grid1 use the same basis).

 %---------------------------------------------------------------------
 % Set default data
 d = grid1{1}.d;
 xi_default = 1.0/(grid1{1}.d+1); % barycenter
 x_default = 0;
 %---------------------------------------------------------------------

 %---------------------------------------------------------------------
 % Interpolate at the quad. points
 igi = 1;
 iei = 1;
 np   = size(xxx,2);
 nvars = size(uuu1{1},1); % number of variables
 uue  = zeros(base1.pk,nvars);
 xi   = zeros(size(xxx,1),1);
 PHI  = zeros(base1.pk,1);
 ui   = zeros(nvars,np);
 tic
 for ip=1:np
   if(mod(ie,floor(np/100))==0)
     toc
     disp(['Element ',num2str(ie),' of ',num2str(np)]);
     tic
   end
   [igi,iei,xii] = locate_point(xxx(:,ip),grid1,ddc_grid1,[igi,iei],1);
   if(iei<0)
     warning(['Unable to locate point ',num2str(xg(:,l)'), ...
           ' in element ',num2str(ie),'; using default value.']);
     xi(:) = xi_default;
     uue(:,:) = x_default;
     igi = 1; iei = 1;
   else
     for id=1:nvars
       uue(:,id) = uuu1{igi}(id,:,iei);
     end
     xi(:) = xii(2:d+1);
   end
   % interpolate
   for i=1:base1.pk
     PHI(i) = ev_pol(base1.p_s{i},xi);
   end
   for id=1:nvars
     ui(id,ie) = sum( uue(:,id).*PHI );
   end
 end
 %---------------------------------------------------------------------

return

